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We study the interactions between two atomic species in a binary Bose-Einstein condensate to 
revisit the conditions for miscibility, oscillatory dynamics between the species, steady state solutions 
and their stability. By employing a variational approach for a quasi one-dimensional, two-atomic 
species, condensate we obtain equations of motion for the parameters of each species: amplitude, 
width, position and phase. A further simplification leads to a reduction of the dynamics into a simple 
classical Newtonian system where components oscillate in an effective potential with a frequency 
that depends on the harmonic trap strength and the interspecies coupling parameter. We develop 
explicit conditions for miscibility that can be interpreted as a phase diagram that depends on the 
harmonic trap's strength and the interspecies species coupling parameter. We numerically illustrate 
the bifurcation scenario whereby non-topological, phase-separated states of increasing complexity 
emerge out of a symmetric state, as the interspecies coupling is increased. The symmetry-breaking 
dynamical evolution of some of these states is numerically monitored and the associated asymmetric 
states are also explored. 



I. INTRODUCTION 

Over the past 15 years, the study of Bose-Einstein con- 
densates (BECs) has gained a tremendous momentum, 
stemming from an intense and wide variety of theoreti- 
cal, as well as experimental studies that have now been 
summarized in a number of books [1, 2]. One of the par- 
ticularly intriguing aspects of the system is its effective 
nonlinearity stemming from a mean-field representation 
of the inter-atomic interactions. This, in turn, has led 
to a wide range of developments in the area of nonlin- 
ear matter waves in BECs [3] and the drawing of natural 
parallels between this field and that of nonlinear optics, 
where similar nonlinear Schrodinger (NLS) types of mod- 
els arise [4]. 

One of the particularly interesting aspects of investi- 
gation of BECs within the realm of NLS (often referred 
to in the BEC context as Gross-Pitaevskii) equations is 
based on the examination of multi-component systems. 
Starting from the early work on ground-state solutions 
d, and small-amplitude excitations [7|] of the order pa- 
rameters, numerous investigations have been focused on 
the study of two hyperfine states or two different atomic 
species that can be condensed and confined concurrently. 
More specifically, a few among the numerous topics in- 
vestigated involved the structure and dynamics of binary 
BECs [8, 9, 10], the formation of domain walls between 
immiscible species , bound states of dark-dark [l2[ , 

or dark-bright [l^, or coupled- vortex [15], or even 
spatially periodic states [l^. The early experimental ef- 
forts produced such binary mixtures of different hyper- 
fine states of ^'^Rb pJi] and of ^^Na [19], but also of mixed 



*URL: 
tURL: 



http : //rohan . sdsu . edu/^rcarrete/ 1 



http : //nlds . sdsu . edu| ~ 



condensates [18|. Efforts were later made to create two- 
component BECs with different atomic species, such as 
'^^K-^'^Rb [20] and ^Li-^^^Cs [21], among others. Re- 
cently the interest in multi-component BECs has been 
renewed by more detailed and more controlled experi- 
mental results illustrating the effects of phase separation 
[22I . I23I . [24| . which have, in turn motivated further the- 
oretical studies in the subject [25|, [2^; see also [27] for a 
recent review. 

Although in the present work, we will focus on two- 
component condensates, it is relevant to note in passing 
the increasingly growing interest in three-component, so- 
called, spinor condensates [28]. Among the numerous 
themes of investigation within the latter context we men- 
tion spin domains [29], polarized states [30], spin textures 
[sij , as well as multi-component (vectorial) solitons of 
bright J32, [H, m, [al, dark ^], gap [s^, and bright- 
dark [3q] types. 

Our aim in the present manuscript is to revisit the 
theme of binary condensates in quasi-one-dimensional 
BEC settings, in an attempt to offer additional both an- 
alytical and numerical insights on the phenomenology of 
phase separation. Our manuscript is organized as fol- 
lows. In Sec. ini a Gaussian trial function is used in 
a variational approach to obtain six first-order ordinary 
differential equations (ODEs) for the time evolution of 
the parameters of the two-component ansatz: position, 
amplitude, width, phase, wave number and chirp. In 
Sec. nil Al the fixed points of the system of ODEs are 
obtained to yield the equilibrium position, amplitude, 
and width of the two species. Bifurcation diagrams of 
the ansatz' parameters are produced as a function of the 
interspecies coupling strength. In Sec. IIIIB[ phase dia- 
grams are produced and analytical conditions for mis- 
cibility are obtained relating the interspecies coupling 
with the system's chemical potential and parabolic trap 
strength. The dynamics of this system is compared to 
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results obtained by numerical integration of the Gross- 
Pitaevskii equation in Sec. IIII Ci In Sec. IIII Dl the sys- 
tem of ODEs is further reduced, upon suitable approx- 
imations, to a classical Newtonian system; the latter is 
simple to analyze and instructive with respect to the in- 
terpretation of the fundamental interactions driving the 
system's dynamics. In Sec. IIV( we numerically analyze 
in a systematic way the existence and stability of higher 
excited phase-separated states as a function of the in- 
terspecies interaction. The dynamical instability evolu- 
tion of the latter class of states is examined numerically 
in Sec. [Vl Motivated by the numerical experiments of 
Sec. |Vl in Sec. I VII we study the existence of asymmetric 
states when the chemical potentials of the two compo- 
nents differ. Finally, in Sec. I VIII we summarize our re- 
sults and present some interesting directions for future 
investigation. 

II. VARIATIONAL MODEL 
A. Coupled equations 



In a highly anisotropic trap, where the frequency of 
the longitudinal component of the trap is much smaller 
than the transverse components cOx ^ oOy = a;^, an effec- 
tive one-dimensional (ID) system of partial differential 
equations (PDEs) can be obtained ^] 

dui / 1 9^ 9 , ,2 , ,2\ 

^^ = [-,9^ + Y^' + \-^f+9\u.f)-u (2) 

where time, space, and wave function have been rescaled 
to reduce the system's parameters to just two (Vt and g). 
In this ID reduction, the chemical potential corresponds 
to jiiD-, i^j Qjj / jiiD^ X x^mjiiD/h is the longi- 

tudinal distance from the center of the trap, t tfim/h 
is the rescaled time, g = ^12/^11 = ^21/^22 is the rescaled 
species interaction term. In the literature, the condition 
of miscibility A = (012021 - ^ii^22)/^ii = - 1 is of- 
ten used see e.g. [sl, |9|, [lQ|- The rescaled harmonic trap 
frequency is given by = huJx//^iD- 



The Gross-Pitaevskii (GP) equation, which is a variant 
of the NLS equation accounting for the potential confin- 
ing the atomic species, governs the dynamics of bosonic 
particles near absolute zero temperatures P, Q . In the 
context considered herein (related to the case of ^^Rb 
which is common in relevant experiments [13, [Hl)? 
two hyperfine states of the same atom are described by 
a set of coupled GP equations [9] 



ih 



dt 



- — + F2 + 922 1^2^ + 921 \^lf ] V^2, 



Vj = -{u;%x^^u;]y^u;%z^)j = l,2. (1) 



Here, 



9jk 



Aiih^ajk/m are the self coupling interaction 



parameters of the first species for j = 1, /c = 1, for the 
second species j = 2, = 2 and for j = 1, /c = 2 or j = 2, 
/c = 1 are the cross species coupling parameters. These 
interaction strengths depend on the scattering lengths 
between same (an and a22) and different species (ai2 = 
a2i). The external harmonic trapping potential for each 
atomic species, Vj = Vj{T) depends on the radial distance 
r from the center of the trap, while the atomic density 
is given by the square modulus of i/jj = il)j{Y^t). The 
atomic mass is denoted by m. 

We will make the customary assumption that the ex- 
ternal trap's effect on each species is the same, making 
Vi{v) = V2{t). More importantly, in the interest of an- 
alytical tract ability of our results, we will also assume 
that the self interactions for each species is the same, 
9ii = 922- Both assumptions are very good approxima- 
tions of the physical reality, although they are not exact; 
see e.g., the relevant discussion of [23]. 



B. Ansatz and Euler-Lagrange Equations 

To develop a variational model, we substitute in the 
rescaled two-component Lagrangian 



where 



f 

J — ( 



(Li + 1/2 + 1/12 + L21) dx, 



duj 



dx 



V{x)\u, 



1 



L21 = ^9\uif \u2f , 



the Gaussian ansatz of the form 



iiC-\-Dx-\-Ex^) 



Ui{x^t) = Ae e 
U2{x,t) = Ae-^^^e^(^-^^+^^ 



(4) 
(5) 



Starred variables indicate complex conjugation. The pa- 
rameters A, 5, C, and W are assumed to be time 
dependent and they represent the amplitude, position, 
phase, wave number, chirp, and width of the Gaussian 
ansatz, respectively. For large 1^, the steady state solu- 
tion very closely resembles two Gaussian functions sep- 
arated by a distance 25, with constant rotation of the 
phase C = /it, where /i is the chemical potential (with- 
out loss of generality we take /i = 1). Upon interaction 
between the species and with the harmonic trap, accel- 
eration induces an inhomogeneity in the carrier wave, 
known as chirp, that accordingly affects the phase of each 
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species. It is important to note here that the ansatz is 
invariant upon transposition of space and atomic species 
component x —x and ui U2. This ahows us to 
reduce the number of parameters in the system to six 
instead of having twelve, which would take into account 
independent variation of the parameters in ui and U2. 
However, this simplification comes at a certain cost as, 
in particular, it is not possible to monitor asymmetric 
(between the two components) states within this ansatz; 
the latter type of states will be partially explored within 
the dynamics of the species in Sec. IVl 



When the Lagrangian is evaluated for the proposed 
ansatz, a spatially-averaged effective Lagrangian is ob- 
tained 



rameters of our ansatz: 



1 ^^W^ 



at at at 



(6) 



and the equations of motion for the parameters are ob- 
tained through the corresponding Euler-Lagrange equa- 
tions 



dA 

~dt 
dB 

~dt 

dC 

~dt 



dD_ 



dE 



dW 

~dr 



-AE, 



D + 2BE, 

12 



(8) 
(9) 



(25^ - 5W^) 



V2A-9 



e"iW^ (85^ + 2B^W^ + 5W^), (10) 



V2A^Bg _ . 
B 

2W^ ' 



V2A^9 

V2A^ 

2EW. 



1 



2 

2E^ 



2 ' 



(11) 

(12) 
(13) 



As described below, these ODEs reflect fairly accu- 
rately the principal dynamical features of the system. 
In particular, they capture the oscillations of the two 
species past each other when the equilibrium separation 
B is zero and the oscillations about their corresponding 
phase-separated equilibrium position when the two com- 
ponents are phase-separated. 



III. STEADY STATE SOLUTIONS 



dL d r dL \ _^ 
dpj dt \dpj J 



where the parameter pj, j = 1,2,..., 6 represents the 
parameters in the ansatz A, B, C, E^ W and pj = 
dpj/dt. The Euler-Lagrange equations for A, 5, and W 
evaluate to dL/dA = 0, dL/dB = and dL/dW = 
since the second term in the Euler-Lagrange equation 
([7j) is zero for these. This results in equations that in- 
volve linear combinations of dC/dt, dD/dt and dE/dt 
that when solved, give Eqs. ^ and (O. The 

rest of the Euler-Lagrange equations give equations that 
involve linear combinations of dA/dt, dB/dt and dW/dt 
that can be solved to give Eqs. (j8|), © and ([T3|). The 
following equations are the result of solving the Euler- 
Lagrange equations for the time derivatives of all the pa- 



A. Phase Bifurcations 

The steady state of Eqs. (|8|)- (p!3|) is obtained by setting 
dA/dt dB/dt dC/dt dE/dt dD/dt dE/dt = 
dW/dt = witch leads immediately to the steady state 
solution E^ = = and = ji. When the equilibrium 
separation between the species is zero, the equilibrium 
amplitude and width reduce to 

5* = 0, (14) 
2^2 f 8/i - ^Ib^t^ + AiiA 

= ^ 

When the equilibrium separation is nonzero, the result- 
ing equilibrium amplitude and width are given by the 




FIG. 1: (Color online) (a) Steady state solution, ui and U2, for 
the mixed state when B = 0. (b) Steady state solution for the 
separated state when B ^ 0. The solid (blue) line represents 
the steady state of the GP equations and the dashed (red) 
line is the steady state solution of the system of ODEs. Here, 
VL 0.6, /i = 1 ^ = 1 for (a) and ^ = 20 for (b). 



transcendental relations 



Bj 



2Wi 



0, (17) 
0, (18) 



mi 



\n''{W^ + 2Bl) = 0. (19) 



Figure [T] shows that the steady state solution of the 
full GP model closely matches the steady state from the 
ODEs, indicating that the ansatz successfully captures 
the relevant PDE behavior. For large values of 1^, the 
steady state solution of the GP deviates from the Gaus- 
sian shape, resembling an inverted parabola, often re- 
ferred to as the Thomas- Fermi approximation [l|, Q . 

Detailed bifurcation diagrams for the amplitude, po- 
sition, and width of each species can be obtained by 
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FIG. 2: (Color online) Equilibrium (a) separation, (b) am- 
plitude, and (c) width for the two condensed species, as a 
function of the interspecies coupling strength g for Q = 0.6 
and \i — \. The (red) thin solid (stable) and dashed (un- 
stable) lines correspond to the steady states for the reduced 
ODE model [Eqs. ©-{[IS))], while the (blue) thick solid (sta- 
ble) and dots (unstable) correspond to the steady state from 
the full GP model [Eqs. (E))-©]. 



solving Eqs. ([I71)-([T9]) for A, B, and W for the phase 
separated state and Eqs. (p!4|) - (p!6|) for the mixed state. 
The steady state solution reveals a pitchfork bifurcation 
for the position of each condensate as the interspecies 
coupling strength is increased as can be seen in Fig. [2^. 
Interestingly, the steady state of the GP equations pro- 
duces a supercritical pitchfork bifurcation at point D of 
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Fig. O More specifically, for small values of ^, a sta- 
ble mixed phase can be identified; as g is increased, the 
mixed phase becomes unstable and the phase separated 
state becomes stable. On the other hand, the system of 
ODEs also predicts a pitchfork bifurcation, however the 
approximate nature of the ansatz results in the identifi- 
cation of the bifurcation as a subcritical one (point A) oc- 
curring in the vicinity of a symmetric pair of saddle node 
bifurcations (points B). Both bifurcation diagrams agree 
with each other away from the transition between phases. 
In the vicinity of the transition point, clearly, the nature 
of the ansatz is insufficient to capture the fine details of 
the bifurcation structure (thus inaccurately suggesting a 
phenomenology involving bistability, and hysteretic be- 
havior of the system). It should be noted that the phase 
separation from the GP model (point D) lies near the 
saddle node bifurcation point (point B) and the subcriti- 
cal bifurcation point (point A) from the system of ODEs. 
At small values of the harmonic trap strength, the true 
point of phase transition is closer to the subcritical bi- 
furcation point and at larger values of the harmonic trap 
strength, the true point of phase transition lies closer to 
the saddle node bifurcation point. 



B. Phase Separation 
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FIG. 3: (Color online) The (blue) solid line represents the 
boundary of zero species separation from the GP equations 
(point D in Fig. [2]), where states that lie to the right of this 
line are mixed and values to the left are phase-separated. The 
(red) dashed line represents the boundary of zero species sep- 
aration for the system of ODEs given by Eq. (|2Q|) (point A in 
Fig. [2]). The (black) dotted line shows the saddle- node point 
(point B in Fig. [2]) obtained from the system of ODEs. 



Using the zero separation amplitude and width, point 
A in Fig. [21 an expression for the onset of phase separa- 
tion can be approximated in terms of the critical inter- 
species interaction 

6, + 3vWTV^ (20) 
26/i - 7v/T5lFTV 

Despite the deviation of point A in Fig. [21 from the rel- 
evant point D of the corresponding PDE, the analytical 
expression offers valuable insight on the dependence of 
the critical interspecies interaction for phase separation 
on parameters of the trap (in particular, its frequency) 
and those of the condensate (in particular, its chemical 
potential); see also Fig. [3] for a detailed comparison of 
the ODE and PDE bifurcation points. More specifically, 
the equation predicts that when the harmonic trap's fre- 
quency approaches zero, ^ 0, phase separation occurs 
when ^cr 1 coinciding with the well known miscibility 
condition A = - 1 = 0. As ^ l^cr = 4V2/i/7 ^ 0.8/i 
phase separation occurs at ^cr oo. This behavior can 
be qualitatively understood since tighter (larger traps 
tend to "squeeze" both components together, thus frus- 
trating the system's tendency towards phase separation. 




FIG. 4: Oscillations for a mixed state from direct numerical 
integration of the GP model {|2])-(|3]) for = 1 (^ < ^cr), M = 1 
and Q = 0.6. Lighter gray corresponds to one the density of 
one component and darker gray to the other component. 



The prediction of the system of ODEs for the loca- 
tion of the bifurcation point agrees well with the results 
from the GP model for small values of the harmonic trap 
strength. Recall that the phase separation for the ODE 
model is located at the subcritical pitchfork bifurcation 
point A in Fig. [21 However, as can be observed from 



Fig. [31 a better approximation for the full system's phase 
separation (point D), in the case of large trapping fre- 
quencies, can be given by the saddle-node bifurcation 
point B for the ODE reduced system for Q values close 

to Qnr- 



6 




FIG. 5: (Color online) (a) Amplitude, (b) position, (c) phase, 
(d) velocity, (e) chirp, and (f) width corresponding to the 
mixed oscillating state show in Fig. U] Solid lines represent 
results from direct numerical simulations of the GP equation 
while dashed lines depict the results for the ODE reduction 

©-(Ull). 




FIG. 6: Same as in Fig. |4]for the time evolution of oscillating 
species about their phase separated configuration when g = 
20 (g > gcr). 




FIG. 7: (Color online) Same as in Fig.[5]for the time evolution 
of oscillating species about their phase separated configura- 
tion as shown in Fig. [6] 



C. Dynamics of the Reduced System 

For relatively small values of the interspecies coupling 
g (i.e., to the left of point D in Fig. [2]) the two species do 
not separate and thus, when given opposite direction ve- 
locities from the mixed state, they will oscillate through 
each other as depicted in Fig. [H This case is analyzed in 
Fig. [5] where the two components oscillate through each 
other about their common equilibrium separation of zero 
and the prediction from the system of ODEs (red dashed 
lines) agrees very well with the direct integration of the 
GP equations (blue solid lines). Because of conservation 
of mass, the amplitude and width oscillate out of phase: 
the amplitude is maximized and width is minimized when 
the acceleration of the species is maximized; on the other 
hand, the width is maximized and the amplitude is min- 
imized when the velocity is maximized. The velocity of 
each species [see Eq. (|9|)] has two components, one that 
depends on the wave number D and another that depends 
on the product of chirp and position 2EB. If there is no 
chirp, the wave number is the velocity of the condensate 
species. In Fig. [St and[7t, a factor of /it has been added 
to show the deviation of the phase from the steady state 
value. 

On the other hand, if we assume a well-separated state 
as the PDEs initial condition, (right of point D in Fig. [2]), 
it is possible that the two components entertain oscil- 



FIG. 8: (Color online) The effective potential, [/gff , for several 
values of the species coupling parameter: g = 0, g = 2.6, 
^ = 6, and g — 20. Here, Q — 0.6 and /x = 1. 



lations about their phase-separated steady state as de- 
picted in Fig. [6l As illustrated in Fig. [71 the two phase 
separated components collide against each other as they 
oscillate (but do not go through each other) about a 
nonzero position. The prediction from the system of or- 
dinary differential equations once again agrees very well 
with the numerical integration of the GP equation for 
the position of the two species, but differs somewhat for 
the other parameters of motion. This occurs because the 
time evolution of the solution of Eqs. ([2]) and (|3]) in this 
case, due to the oscillation and interaction, deviates from 
a Gaussian waveform and the corresponding variational 
prediction begins to lose accuracy. 



D. Newtonian Reduction 

To develop a more tractable model for the dynamics, a 
classical Newtonian system for the motion of the center 
of each species is desirable. Taking a time derivative of 
Eq. dH and substitution of Eqs. jH), ([II]), and ^ yields 



= - 



V2A: 



9 



(21) 



where we simplified the dynamics by replacing the time 
dependent variations of A and W by their respective equi- 
librium values and W^. This simplification is justified 
by the fact that the oscillations in A and W are relatively 
weak as it can be observed from Figs. [5^, [Sf, [7^, and 
[Tf. These phase separated oscillations contain two funda- 
mental physical features: the external trapping potential 
and an exponential repulsive interaction that depends 
on the cross species coupling parameter g. Integrating 
Eq. ([2T]) with respect to B yields a Newtonian equation 
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FIG. 9: (Color online) (a) Amplitude of oscillation and (b) 
period of oscillation, T, as a function of increasing initial ve- 
locity, t'o. The (blue) solid line represents results from the 
coupled GP equations, the (red) dashed line represents re- 
sults from the system of ODEs, and the black dotted line 
represents the results from the Newtonian reduction. Here, 
g — \^ [I — \ and Q = 0.6. 



of motion under the effective potential 



t/eff = \b^ 



V2AI 



(22) 



This reduced dynamics gives an effective double well po- 
tential for a fully phase separated state {g > gcr) and 
a nonlinear single well potential for the mixed states 
{g < 9cr)' It is important to note that in the expres- 
sion of Eq. ([22|). A* and B^ vary as a function of g. This 
dependence has been incorporated in Fig. [8] which shows 
the effective potential for increasing values of ^, where it 
transitions from a single well potential to a double well 
potential. For Q = 0.6 and /i = 1, the reduced Newtonian 
model predicts that at ^ = 2.6 the two species phase sep- 
arate and the potential becomes very flat yielding very 
long periods of oscillation. It is remarkable to point out 
that a similar picture has been drawn for the interaction 
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FIG. 10: (Color online) (a) The period of oscillation, T, is 
shown as a function of g for an initial velocity of vo — 0.01. 
(b) The period of oscillation, as a function of g for an initial 
velocity of vq — 0.1. The (blue) solid line represents results 
from the coupled GP equations, the (red) dashed line rep- 
resents results from system of ODEs, while the black dots 
represent the results from the evaluation of the eigenvalues 
of the Jacobian of the system of ODEs at equilibrium. Here, 
Q — 0.6 and = 1. 



of two particle-like excited states (i.e., dark solitary mat- 
ter waves) within the same species and has been found to 
be extremely successful in comparisons with experimen- 
tal results [39]. In that case, as well, the fundamental 
characteristics were the parabolic trapping and the ex- 
ponential tail-tail interaction between the waves. 

Figure [9] depicts both the (a) amplitude and (b) period 
of oscillations predicted by the Newtonian reduction in 
comparison with the corresponding PDE findings. It is 
seen that the GP equation and the ODE system agree 
very well for a wide range of values of initial veloci- 
ties. This figure also shows the nonlinearity of oscillations 
where small amplitude oscillations have a period of 18.3 
and large oscillations yield a period oi T ^ 10.5, which 
corresponds to the harmonic trap's period T 2T:/Vt. 



Figure [TOk shows that the period of oscillation from the 
system of ODEs matches that of the GP model for small 
values of g. As g ^ 2.4 (i.e., when approaching phase 
separation of the GP equations [point D in Fig. [2]), the 
period changes substantially. Then as g 3.1, the sys- 
tem of ODEs begins to phase-separate. Since these are 
small oscillations, the eigenvalues from the Jacobian of 
the system of ODEs at equilibrium matches very well 
the oscillations of the system of ODEs. Figure [TOb shows 
similar results to Fig. [TOk but for larger oscillations. We 
can see that for larger oscillations the ODEs' period more 
closely matches that of the GP for smaller ^, while for 
larger values of the interspecies strength, the deviation 
becomes more significant. Furthermore, for larger val- 
ues of g, the period obtained from the eigenvalues also 
deviates from the period from the ODEs. For large oscil- 
lations, the two species do not effectively interact (since 
they go through each other too rapidly to "feel" each 
other) and the period is roughly independent of g as pre- 
dicted by Eq. ([2T]) when \B\ <C 1, yielding simple har- 
monic oscillations. In this case, the GP's period is close 
to that predicted by the ODEs. 

IV. EXCITED STATES 

For sufficiently weak traps Q < 0.5, as the interspecies 
coupling strength is increased for ^ > 1, new, high order 
mixed states emerge. These higher excited states cor- 
respond to alternating bands dominated successively by 
each of two species. As ft is decreased or alternatively 
g is increased, the number of alternating bands increases 
within the solution profile and the population imbalance 
within each band less pronounced. Each solution branch 
is found by using parameter continuation on the param- 
eter g using a Newton fixed point iteration to find the 
stationary state. The stability for each computed pro- 
file was obtained by computing the eigenvalues of the 
linearized dynamics at the fixed points using standard 
techniques |3i]. We now describe the series of bifurca- 
tions that occur as the interspecies coupling g increases 
as depicted in Fig. [TTJ 

• For g < 1.08, the only state that exists is the mixed 
state and it is stable. This threshold is close to the 
traditional miscibility condition A = — 1 = 0. 

• At ^ = 1.08 (point A), the mixed states becomes 
unstable, through the previously described super- 
critical pitchfork bifurcation, leading to the emer- 
gence of a stable phase-separated state where one 
component moves to the left and the other one 
moves to the right. 

• At ^ = 1.24 (point B), a second supercritical pitch- 
fork occurs, rendering the symmetric (across com- 
ponents) solution more unstable and giving rise to 
a state where one species has a single hump and the 
other one has a double hump. We call this state a 
1-2 hump configuration. 
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FIG. 11: (Color online) Degree of phase separation Ah = ui(xc) — U2{xc), where Xc is the location of the closest density 
maximum (irrespective of component) to the trap center, as a function of the interspecies coupling parameter for various 
excited states for g G [0.5, 4.5] and Q = 0.2. Stable (unstable) solutions branches are depicted with black solid (dashed) lines. 
Typical solutions for each branch are depicted in the surrounding insets. 




g 



FIG. 12: (Color online) Real part of largest eigenvalue of the 
mixed state as a function of the interspecies coupling param- 
eter for the various excited states depicted in Fig. [TT] with 
g G [0.5,4.5], and Q = 0.2. 



• At ^ = 1.63 (point C), the third supercritical pitch- 
fork bifurcation of the series arises, leading this 
time to a situation where one species moves to the 
left and the other one moves to the right, both 



forming a double hump (a 2-2 hump configuration). 

• At ^ = 2.22 (point D), a double hump with a triple 
hump state arises (a 2-3 hump configuration). 

• At ^ = 3.8 (point E), a triple hump with a triple 
hump state forms (a 3-3 hump configuration). 

As the condensate becomes wider, more and more bands 
appear, each band having a width comparable to the 
healing length of the condensate. This picture seems to 
be natural from the point of view of [25[ , where the emer- 
gence of these higher excited states could be interpreted 
as a manifestation of an effective modulational instabil- 
ity. Within the effectively finite region determined by the 
confining potential, as g is increased, higher "modulation 
wavenumbers" become unstable, leading to the "quan- 
tized" (associated with the quantization of wavenum- 
bers in the effectively finite box) cascade of supercritical 
pitchfork bifurcations and associated further destabiliza- 
tions of the symmetric state. The degree of phase sep- 
aration in Fig. [11] between these bands is computed as 
Ah = ui{xc) — U2{xc)^ where Xc is the location of the 
closest density maximum, irrespective of component, to 
the trap center. 

The emergence of high order states can be inferred 
by observing the real part of the eigenvalues as the in- 
terspecies coupling is varied. The eigenvalues collide 
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with the imaginary axis as new states emerge. Eigen- 
value analysis shows that all excited states are unstable 
with the exception of the single-hump phase-separated 
state (which results from the first pitchfork bifurcation, 
namely the 1-1 hump state). These excited states are, 
however, less unstable than the mixed state from which 
they arise. The latter, as shown in Fig. [121 becomes pro- 
gressively more unstable, as expected, as further multi- 
hump branches arise. Each of the bifurcation points, for 
which the same designation as in Fig. [11] is used, corre- 
sponds to a further pair of real eigenvalues appearing for 
the mixed branch. 



o 


1-2 hump 


E '\, 


2-2 hump 




2-3 hump 




3-3 hump 


D 












O N 

\ 'n, 

\ ^ . 




\ ' ^ _ 





^0 5 10 15 20 25 

g 



FIG. 13: (Color online) Real part of the largest eigenvalues 
of the principal phase-separated states as a function of the 
interspecies coupling parameter, ^, for — 0.2. 




As indicated by Fig.[TTl the 1-1 hump phase separated 
state is stable even for large values of ^. In fact, for 

<^ 1 and \^ the two components repel each other 
so strongly that the center of magnetic trap becomes 
a domain wall, where each species abruptly transitions 
from near zero atomic density to maximum atomic den- 
sity. Apparently, this 1-1 hump phase-separated state 
is the only stable state of the system after the mixed 
state loses its stability past the bifurcation point A (see 
Fig. [TT]) . Nonetheless, as depicted in Fig. [131 the addi- 
tional multi-humped excited states are successively cre- 
ated as g is increased (bifurcation points B, C, D, and 
E in Fig. [11]) and are progressively less unstable as g 
increases. Therefore, for sufficiently large ^, the instabil- 
ity of higher excited multi-humped states might be weak 
enough for these states to be observable within experi- 
mentally accessible times. It is interesting to note that 
since these excited multi-humped states emanate from 
the (already unstable) mixed state, they feature a rela- 
tively strong instability close to their bifurcation point. 



FIG. 14: (Color online) Dynamics of the mixed state for g — ^ 
and Q = 0.2. (a) Top and bottom subpanels correspond to 
the evolution of the densities for the first and second compo- 
nents, respectively, after applying an initial spatially random 
perturbation of size 1 x 10~^. Panels (b) and (c) depict snap- 
shots of the initial density and at the times indicated in panel 
(a) by the white vertical dashed lines, (d) Long term dy- 
namics showing that the mixed state eventually approaches a 
separated state, (e) Snapshots of the densities at t = and 
t — 10, 000 (see white vertical dashed line in panel (d)). 



V. DYNAMICS OF UNSTABLE STATES 

In this section we present the dynamics of the differ- 
ent unstable states that were identified above. We start 
by analyzing the destabilization of the mixed state for 
a value of g to the right of the bifurcation point A (see 
Fig. [TT]) . In fact, for all the dynamical destabilization 
results presented in this section we chose a value of ^ = 5 
and 1] = 0.2 that is to the right of the bifurcation point 
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FIG. 16: (Color online) Dynamics of the unstable 2-2 hump 
state for ^ = 5 and Q = 0.2 as in Fig. 1141 The eventual result 
of the instability is the formation of a robust and dynamically 
stable 1-1 hump state. 



FIG. 15: (Color online) Dynamics of the unstable 1-2 hump 
state for ^ = 5 and Q = 0.2. (a) Same as in Fig. (Mja. (b) 
Snapshots of the densities at t = and t = 758 (see white 
vertical dashed line in panel (a)), (c) Growth of the norm of 
the perturbation vs. time (in semi- log plot). Circles (crosses) 
correspond to the perturbation of the first (second) compo- 
nent and the red solid line is the instability growth obtained 
from the eigenvalue computation (max(3f^(A))) = 0.122). 



E (see Fig. [TT]) . Therefore, in this regime we have co- 
existence of several unstable multi-hump solutions and 
the stable phase-separated 1-1 hump state. Other pa- 
rameter regions (not shown here) gave qualitatively the 
same results. 

In Fig. [m we show the destabilization of the mixed 
state. As it can be observed from the figure, the mixed 
state suffers an initial modulational instability that be- 
comes apparent for t > 35 that seeds a highly perturbed 
2-3 hump solution [see panel (b)]. Since this new solu- 
tion is also unstable for the chosen parameter values, it 
is rapidly converted into a relatively long lived 1-2 hump 
state [see panel (c)]. Nonetheless, as it is clear from the 
long term dynamics presented in panel (d), the 1-2 hump 
state, being unstable, eventually "decays" to the sepa- 
rated (1-1 hump) state. We use here the term "decaying" 
in quotes since our system is conservative (no dissipation) 
and thus there is no real decay. 

In Fig. [15] we show the destabilization of the unstable 
1-2 hump state. As it is obvious from the figure, al- 
though the instability eigenvalue is sizeable [A^ = 0.122, 
see panel (c)]), the initially weak perturbation does not 
lead to the destruction of the 1-2 hump state. Instead, 
this unstable state just momentarily "jumps" to the left 
(or to the right, results not shown here), see panel (b). 
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FIG. 17: (Color online) Dynamics of the unstable 2-3 hump 
state for g — 5 and Q — 0.2 as in Fig. 1141 The initial in- 
stability leads to a transient 1-2 state that eventually ap- 
proaches the separated state for extremely long propagation 
times (t > 100,000). 
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FIG. 18: (Color online) Dynamics of the unstable 3-3 hump 
state for g — 5 and Q — 0.2 as in Fig. 1171 The initial insta- 
bility dynamics reshapes the waveform into a 2-2 state, then 
to long-lived transient a 1-2 state that eventually approaches 
the 1-1 state for t > 40, 000. 



for a short period of time (becoming slightly asymmetric) 
and then comes back close to the 1-2 hump (symmetric) 
steady state configuration. This indicates that, for this 
parameter combination, the 1-2 hump state is a saddle 
fixed point. Thus, the orbit remains close to the steady 

1- 2 hump state and it is eventually "kicked out" along the 
unstable manifold. Then, it performs a quick excursion 
and "comes back" through the stable manifold. 

In Fig. [16] we depict the destabilization dynamics of the 

2- 2 hump state. As it is evident from the figure, when 
compared to Fig. [151 the destabilization happens earlier 
{t ~ 200) since the 2-2 hump state is more unstable than 
the 1-2 hump state (see Fig. [13]). As the dynamics of 
the perturbed 2-2 hump state evolves, it progressively 
favors phase separation until eventually the system rear- 
ranges itself into a 1-1 hump state at about t ^ 4000. 
Since the 1-1 hump is stable, this resulting configuration 
is sustained thereafter. 

A similar scenario (fast initial destabilization, slow 
transient stage and eventual settling into the stable 



phase-separated state) is observed when following the 
dynamical destabilization of the 2-3 hump state. As it 
can be observed from Fig. [171 the relevant configuration 
destabilizes around t ~ 90 (earlier than its 1-2 and 2- 
2 hump counterparts given its larger instability growth 
rate [see Fig(T3]). The dynamics goes through a transient 
1-2 state and eventually settles into a highly perturbed 
separated state that is preserved thereafter as was the 
case also for the mixed state dynamics (see Fig. [14]) and 
the 2-2 state (see Fig. [T6|) presented above. 

Finally, in Fig. [18] we present the dynamical destabi- 
lization for the 3-3 hump state. Again, this state desta- 
bilizes even faster (t = 40) than the previous waveforms 
because of its higher instability eigenvalue (see Fig. [T3]). 
In this case, the 3-3 hump state first destabilizes into a 
highly perturbed 2-2 state. Since the 2-2 hump state is 
unstable, the ensuing dynamics results into a separated 
1-2 hump state after t > 100 which, in turn, eventually 
settles to a highly perturbed separated state. 
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FIG. 19: (Color online) Dynamics of the real part of the wave 
functions for the 2-3 hump state of Fig. [17] for ^ = 5 and 
Q = 0.2. Panels (a) and (b) correspond, respectively, to times 
centered about t — 50, 000 and t — 100, 000. The blue square 
(green circles) correspond to the first (second) component. 
The continuous lines correspond to the best sinusoidal fit from 
which the local chemical potential can be extracted. 



VI. EXISTENCE AND STABILITY OF 
ASYMMETRIC STATES 

It is interesting to observe in the above computations 
that the transient 1-2 hump state emanating from the 
destabilization of the mixed state (see Figs. [H|) or from 
higher order states (see Figs. [17] and [TH]) has a relatively 
long life span. However, being this state also unstable, 
it eventually tends to a perturbed separated (i.e., 1-1 
hump) state. The process of converting a 1-2 state into 
the separated state involves the effective shift of mass 
from one of the two humps (in the component with two 
humps) to the other one until all the mass is "swallowed 
up" by one hump resulting in a 1-1 hump. For exam- 
ple, as it can be observed in panel (d) of Fig. [ITl at 
t = 50,000 the right hump of the first component has 
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FIG. 20: (Color online) Local chemical potential at the loca- 
tion of maximum density for each component. Blue solid line 
corresponds to the first component while the green dashed line 
corresponds to the second component. The different panels 
correspond to the evolution of the local chemical potentials for 
the steady states: (a) mixed state (cf. Fig. [M), (b) 2-3 hump 
state (cf. Fig.ini, and (c) 3-3 hump states (cf. Fig. fTSj) . 



more mass than its left hump. Since the chemical po- 
tential (i.e., rotation frequency of the wave function in 
the complex plane) is closely related to the mass of the 
condensate, it is possible to follow the local change in 
mass by following the local chemical potential. In Fig. [19] 
we depict the oscillations of the real part of the wave 
functions at the location of the maximal density (see 
squares and circles in the figure). We then fit a sinu- 
soidal curve (see solid lines in the figure) through the 
data to obtain the oscillation frequency. As it is obvious 
from the figure, although both components started with 
the 2-3 hump steady state with the same chemical po- 
tential /ii = /i2 = /i = 1, the two components oscillate at 
different rates. In fact, at t = 50, 000 the local chemical 
potentials for each species is, respectively, jui = 0.776 and 
= 0.734, and at t = 100, 000 we have /ii = 0.688 and 
fi2 = 0.771. We systematically extracted the local chem- 
ical potentials at the location of the maximum density 
for the mixed state, 2-3 hump state and 3-3 hump state 
and depict them in Fig. [20l As it can be observed from 
the figure, both chemical potentials start at /ii = /i2 = 1 
but rapidly drop when the initial instability of the steady 
state develops for the three cases. Then, the chemical po- 
tentials slowly drift until they acquire values very close to 
each other at around the time when the dynamics settles 
to a perturbed phase separated state. 

The fact that during the dynamical evolution of the 
unstable higher order states the local chemical poten- 
tials vary naturally prompts the question of existence and 



stability of steady states with different chemical poten- 
tials between the species. Up to this point all the steady 
states we analyzed supposed that both components had 
the same chemical potential. Namely, the evolution for 
each components is 

U2{x,t) = W2{x)e-'^'\ 

where the steady state profiles are wi^2{x) and /ii = /i2 = 
1. We now relax this symmetric assumption and look for 
steady states with different chemical potentials between 
the components (/ii 7^ /i2)- A full bifurcation diagram as 
a function of /ii, for fixed /i2 = 1, is depicted in Fig. [211 
In this diagram we use in the vertical axis the difference 
in the variance of the steady state distribution between 
the two components. The different bifurcation branches 
correspond to the particular asymmetric states depicted 
in the surrounding insets. This diagram provides an ac- 
count of how the different solutions bifurcate from each 
other. It is very interesting to follow the bifurcation path 
of all the states labeled from A to T during which all the 
higher order states are browsed continuously. Due to 
symmetry, a similar scenario is present when one keeps 
/ii = 1 fixed and varies /i2 (results not shown here). 

The stability in the bifurcation diagram depicted in 
Fig. [21] is denoted, as before, with a solid line for stable 
states and a dashed line for unstable states. As it can be 
noticed from the figures, as it was the case for symmetric 
states, all the asymmetric states are unstable except the 
mixed state. This is due to the fact that the interspecies 
coupling for these diagrams was chosen as ^ = 5, i.e., 
high enough so that the mixed state in unstable. In fact, 
in Fig. [22] we show the largest real part of the eigenvalues 
for the different asymmetric states depicted in the bifur- 
cation diagram in Fig. [211 As it was observed for the 
symmetric case, the higher the order of the asymmetric 
state the more unstable it becomes. It is crucial to note 
that the 1-2 hump state seems to get stabilized for /i2 = 1 
and /ii > 1.12 (see top panel of Fig. [22]) . However, af- 
ter close inspection, see the bottom panel in Fig. [22l the 
real part of the eigenvalue for the 1-2 hump state never 
vanishes but instead becomes extremely small (between 
10~^ and 10~^) until the branch disappears (at jui ^ 1.65 
for /i2 = 1)- This very weak instability explains the fact 
why the transient 1-2 hump state appears to be sustained 
for very long times in the evolutions depicted in Figs. [HI 
[13 and[l8l 



VII. CONCLUSIONS 

In the present work, we have analyzed the emergence 
of non-topological, phase-separated states in the immis- 
cible regime out of mixed ones in the miscible regimes, 
as a natural parameter of the system (namely the inter- 
species interaction strength) was varied. Our analysis 
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FIG. 21: (Color online) Bifurcation diagram of asymmetric states for constant /i2 = 1 as a function of /ii. The vertical axis 
corresponds to the difference of the variance of the steady state configurations for both components for g = 5 and Q = 0.2. 
Stable (unstable) solutions branches are depicted with black solid (dashed) lines. Typical solutions for each branch are depicted 
in the surrounding insets. 




FIG. 22: (Color online) Real part of largest eigenvalue for 
the asymmetric mixed states as a function of /xi with /i2 = 1, 
g = b and Q = 0.2. The bottom panel corresponds to the 
eigenvalue for the 1-2 hump state plotted in logarithmic scale. 



was presented for the case of magnetically trapped two- 
component Bose-Einstein condensates i.e., two incoher- 
ently coupled NLS equations with a parabolic potential. 
Using a variational approach, we obtain a reduced model 
describing the statics, stability and dynamics of each con- 
densate cloud. We are able to elucidate the miscibility 
boundary for the interspecies coupling parameter as the 
strength of the magnetic trap (and/or the chemical po- 
tential of the system) is varied. The approach is also 
capable of accurately capturing the spatial oscillations 
of the clouds about the stable stationary states (both 
for mixed and for phase-separated states). In particular, 
for relatively small interspecies coupling, the two con- 
densed clouds do not phase-separate (mixed state) giv- 
ing rise, if the BEG clouds are initially displaced with 
respect to each other, to oscillations through each other. 
On the other hand, for relatively large interspecies cou- 
pling, the two clouds form a stable phase-separated state 
which can entertain oscillations about the equilibrium 
separation between the components. A further dynam- 
ical reduction allows to understand this behavior more 
intuitively based on an effective potential that undergoes 
a bifurcation from a single well (mixed state) to a dou- 
ble well (phase-separated state) form as the interspecies 
coupling parameter is increased. We also describe the bi- 
furcation scenario of higher order phase-separated states. 
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as the interspecies coupling parameter is increased. We 
observe that several (interwoven between the two compo- 
nents) bands of density modulations progressively arise 
out of the mixed states giving rise to higher excited 
states. Among all the phase-separated states, only the 
first excited one with one hump in each component is 
found to be dynamically stable for all values of the in- 
terspecies interaction strength past its bifurcation point. 
We furthered our analysis by studying the existence and 
stability of asymmetric states for which the chemical po- 
tentials for each species is different. We found, similar 
to the symmetric case of equal chemical potentials, that 
the only stable steady state for high enough interspecies 
coupling is the separated state. Nonetheless, we found 
that in some regimes, the (asymmetric variant of the) 
state with one band in one component surrounded by two 
bands of the other component (the 1-2 hump state) can 
have an extremely weak instability and thus facilitating 
its potential observability in numerical experiments. 

There are numerous directions along which the present 
work can be naturally extended. For example, within the 
one-dimensional context, it is natural to seek to relax 
the simplifying assumption of the intraspecies scattering 
lengths (and by extension the self-interaction coefficients 
of the two components, gn and ^22) being equal. How- 
ever, this extension will unfortunately have to come at 



the expense of an ansatz with different amplitude, width, 
etc. parameters between the components, which will ren- 
der the intuitive and explicit analytical results obtained 
herein much more tedious. 

Another natural extension is to try to generalize the 
ideas presented herein into higher-dimensional or larger 
number of component (i.e., spinor) settings. Especially 
in the former case, more complicated waveforms such 
as crosses and propellers have been predicted in two- 
component condensates in two-dimensions [40] and it 
would be particularly interesting to examine whether 
these, as well as more complicated multi-hump wave- 
forms emerge systematically from the corresponding 
mixed state via two-dimensional generalizations of the 
bifurcations presented herein. Such studies are presently 
in progress and will be presented in future publications. 
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